Use of the reversible jump Markov chain Monte Carlo algorithm to select multiplicative terms in the AMMI-Bayesian model

The model selection stage has become a central theme in applying the additive main effects and multiplicative interaction (AMMI) model to determine the optimal number of bilinear components to be retained to describe the genotype-by-environment interaction (GEI). In the Bayesian context, this problem has been addressed by using information criteria and the Bayes factor. However, these procedures are computationally intensive, making their application unfeasible when the model’s parametric space is large. A Bayesian analysis of the AMMI model was conducted using the Reversible Jump algorithm (RJMCMC) to determine the number of multiplicative terms needed to explain the GEI pattern. Three a priori distributions were assigned for the singular value scale parameter under different justifications, namely: i) the insufficient reason principle (uniform); ii) the invariance principle (Jeffreys’ prior) and iii) the maximum entropy principle. Simulated and real data were used to exemplify the method. An evaluation of the predictive ability of models for simulated data was conducted and indicated that the AMMI analysis, in general, was robust, and models adjusted by the Reversible Jump method were superior to those in which sampling was performed only by the Gibbs sampler. In addition, the RJMCMC showed greater feasibility since the selection and estimation of parameters are carried out concurrently in the same sampling algorithm, being more attractive in terms of computational time. The use of the maximum entropy principle makes the analysis more flexible, avoiding the use of procedures for correcting prior degrees of freedom and obtaining improper posterior marginal distributions.


Introduction
In plant breeding programs, the selection and recommendation of superior genotypes are of considerable importance. The observed performance (phenotype) is determined by the genetic two principal components) representation since the first axes are predicted in a more pronounced manner. A difficulty faced in the use of the Bayes factor or model selection methods that use information criteria is the great demand for time and computational resources. This is due to the fact that all possible models, depending on the dimension, must be adjusted for comparison. From a Bayesian point of view, the uncertainty regarding the best model, which fits a given set of data, can simply be incorporated into the inference problem. Consequently, the model itself is considered as another unknown parameter to be estimated and, with that, each model has its probability distribution. In this context, approaches involving the construction of transdimensional Markov chains, which allow the alteration of the model's dimension during the sampling process, have gained popularity [28].
Among the class of transdimensional methods, there is the Reversible Jump in Markov chain Monte Carlo (RJMCMC), initially proposed by Green [29]. This procedure can be seen as a generalization of the Metropolis-Hastings algorithm when the model's parametric space does not have a fixed dimensionality [30]. This method allows for the flexibility to "jump" between models with parametric spaces of different dimensions, which ensure that the chain is irreducible. Consequently, from the output of a single Markov chain sampler, a complete probabilistic description of the posterior distributions of each model is obtained [28,31,32].
Therefore, the main objective of this study was to adjust the AMMI model under the Bayesian perspective, assuming that the number of principal axes to be retained in the model is a random variable. This assumption was implemented using the RJMCMC methodology for the parameter sampling process. Thus, the joint inference about the number of bilinear terms responsible for the GEI pattern is directly incorporated by the inclusion (or exclusion) of multiplicative terms in the model during the sampling process. Furthermore, specific objectives of the study were as follows: i) to evaluate the sensitivity of the AMMI-Bayesian model to the prior specification for the variance component associated with singular values; ii) use information criteria for model selection, and iii) an evaluation the predictive capacity of the model in a scenario of random loss of genotypes in environment.
To exemplify the use of the method in the evaluation of cultivars we will use simulated and real data. The real data comprise 50 maize genotypes evaluated in 10 locations in southern Brazil. We will address two points that are important in the final stages of plant breeding programs. The first will be to verify the existence of genotypes that have little contribution to the GEI (stable), but that are interesting in terms of yield (broad recommendation). The second will be to identify positive combinations between genotypes and environments, aiming to take advantage of the positive effect of the GEI (recommendation of genotypes to specific environments).

Simulated data
To exemplify the proposed method, a data set with 20 genotypes (G1-G20) and nine environments (E1-E9) was simulated, using a randomized block design, with three replicates. The main effects of genotypes (g) were sampled from a Gaussian distribution with mean equal to zero and variance equal to 12, N(0, 12), and the main effects of the environment sampled from a Gaussian distribution with mean equal to zero and variance equal to one N(0, 1).
The ge ij effects thus simulated make up a GE g×e matrix which is finally corrected for row and column mean so that it corresponds merely to the matrix genotype-by-environment interaction, where g represents the number of rows and e represents the number of columns.
Genotypes belonging to the 1st and 2nd groups were assumed to be unstable genotypes (contribute to interaction). By contrast, the genotypes belonging to the 3rd group were assumed as being stable (they do not have expressive contributions to interaction).
Thus, each value y ijk was obtained by adding the simulated effects plus the error which, in turn, was simulated from a Gaussian distribution N(0,2.5). The simulation was performed using R software [33].

Real data
The real data comprised 50 simple maize hybrids evaluated for grain yield character (t/ha -1 ) during the 2013 and 2014 crop years at 10 locations in the Southern Region of Brazil. The experiment was conducted in an incomplete block design with a variable number of repetitions.
Statistical model. The AMMI model, in vector notation, can be described as: where y (n×1) is the vector of observations with n = ger, in which g corresponds to the number of genotypes, e to the number of environments, and r to the number of repetitions; β (er×1) and g (g×1) are the vectors of block effects within environments and genotype effects, respectively; l k ; a k g�1 ð Þ and γ kðe�1Þ correspond to, respectively, the singular value and the singular vectors (genotypic and environmental) associated with the k-th principal component, with k = 1, . . .,t, in which t is the rank of the GEI matrix, where the singular vectors satisfy the orthonormality constraints and the singular values satisfy the order constraint λ k � λ k+1 � 0. The matrices X 1ðn�erÞ ; X 2ðn�eÞ and, Z (n×g) are designs associated with effects β, γ k and g, respectively. As can be seen, X 1ðn�erÞ and X 2ðn�eÞ are different, where X 2ðn�rÞ is used to distribute the effects of the GE (g×e) interaction and the number of columns is, therefore, equal to the number of environments. The X 1ðn�erÞ matrix, in turn, distributes the effects of blocks within environments, as in Silva et al. [24] and Edwards and Jannink [7] to maintain the identifiability of all parameters and not use marginal conditions where effects sum to zero as presented in Crossa et al. [11]. The vector � (n×1) is composed of random errors with � e N 0; s 2 e I n À � .
If we denote by GE the interaction matrix, then it will have dimensions g × e (g row = genotype and e columns = environment). In the classic AMMI model, the mean of the cell is given by where now α ik and γ jk are coordinates of singular vectors of the k-th principal component. As our model is for parcels, the multiplications involving Z and X 2 in the sum of the model (1) allow distributing and replicating the effects of the GEI in the vector y.
In model (1), the term k specifies the number of axes that are included in the model, which would be from k = 1, where only one axis is being incorporated into the model, up to k = t that comprises the complete model. However, in the approach presented here, the model space is t + 1, such that one can assume an index k = 0, where the corresponding model would be y = X 1 β + Zg + �, that is, the model disregarding the bilinear part.
The likelihood function is given by: Prior distributions for the model parameters are as follows: Hypothesis a: It is assumed that s 2 l k ¼ 1 � 10 8 , which corresponds to assuming a non-informative prior for λ k . This is the hypothesis presented by Oliveira et al. [13], being referred to by Bayesian AMMI (BAMMI).

Hypothesis b:
It is assumed that s 2 l k is a scaled inverted chi-square prior with scale parameter equal to zero, and degrees of freedom equal D ¼ with the restriction of 0 < Δ < 1//2, resulting in a posterior with a shrinking effect for singular values. This is the Bayesian AMMI Shrinkage (BAMMIS) version proposed by Silva et al. [24]: Hypothesis c: It is assumed that s 2 l k is an inverse gamma prior with degrees of freedom equal to one (a = 1) and scale parameter equal to zero (b = 0). This choice is based on the concept of Maximum Entropy (ME), and was implemented in the Bayesian GGE model by Oliveira et al. [27] (more details in S1 Appendix) as follows: As specified, priors are non-informative for β; s 2 e , and genotypic and environmental singular vectors. However, for g, a hierarchical prior at two levels was assigned considering uncertainty about the hyperparameter s 2 g , which results in considering a common population for genotypes. This is similar to what occurs in mixed models when genotype effects are considered random.
The prior distributions, assigned to the singular vectors, are uniform spherical in the corrected subspace due to the orthonormality constraint (the justification for such an assignment is elucidated in Viele and Srinivasan [9]). In this study, there is a special interest in the priors attributed to the singular values hyperparameters. The prior distribution of λ k is a truncated normal for positive values, it is denoted by N + , with zero mean and with variance s 2 l k , which was assumed in the three hypotheses about the prior that was utilized.
Furthermore, the number of principal axes (t) present in the model is a random variable whose a priori uncertainty was established through a truncated Poisson distribution with mean μ. For the μ hyperparameter, a Gamma distribution with scale parameter and degree of freedom equal to one was assigned, as assumed in Balestre and Souza [34].
Posterior complete conditional distributions for model parameters. Combining the likelihood function with the established prior distributions, using Bayes' theorem, we obtain the joint posterior distribution.: The posterior complete conditional distributions are obtained from the joint distribution (3). They are presented, along with the algebraic details in the S2 Appendix.

RJMCMC algorithm decision rule
The objective of using the RJMCMC algorithm was to conduct model selection during the MCMC process. This method allows one to switch between models with different dimensions such that the properties inherent in the MCMC process are not violated. As in this RJMCMC approach, the models are nested, and the decision rule is similar to the proposals of Waagepetersen and Sorensen [35], Bodin and Sambridge [36], and Balestre and Souza [34], where the Jacobian determinant of the transformation in the decision rule is equal to one. Therefore, the decision rule to test the addition of a new main axis to the model is given by: where t corresponds to the current state model and t + 1 to the candidate model with the addition of an axis; p(yt) and p(yt + 1) are the probability functions of data from model t and (t + 1), respectively; p(tμ) and p(t + 1μ) are the priors referring to the respective models; ξ(t,t + 1) and ξ(t + 1,t) correspond to the proposed Hastings correction necessary to guarantee the reversibility condition of the model during the MCMC process, being x t; t þ 1 ð Þ ¼ 1 tþ1 pd and ξ(t + 1,t) = p 0 wherein pa, pd and p 0 correspond to the prior probabilities of adding, deleting, or maintaining, respectively, the number of main axes in the model (pa = pd = p 0 = 1/3).
Simplifying Eq (4), the decision rule for adding a dimension to the model is: If a(t, t + 1) is greater than a random variable u~U(0,1), then the candidate axis is included in the model; otherwise, the model is retained with the current number of main axes present in the model t.
The next step in the algorithm is to carryout the test between the model with dimension t and the model with dimension t − 1. That is, the model with t principal components and the model with reduction of a principal component t − 1, where the decision rule is given by: Simplifying Eq (5), the decision rule to accept or not accept the exclusion of a dimension in the model is: where p(y|t − 1) is the conditional distribution of data for the dimension-reduced model (t − 1) and with ξ(t, t − 1) and ξ(t − 1, t) corresponding to the necessary Hastings corrections to guarantee the reversibility conditions in the sampling, with ξ(t, t − 1) = p 0 and ξ(t − 1, t) = pd/t.

Sampling process
The sampling of the parameters of the models was performed iteratively, combining the Gibbs sampler and the RJMCMC method. The Gibbs algorithm was arranged as follows: Gibbs algorithm. A1-Initial values are assigned to model parameters: A2-The l-th iteration can be obtained as follows:

PLOS ONE
A. Generates β l jg lÀ 1 ; l k lÀ 1 ; s 2 l k � � lÀ 1 ; α k lÀ 1 ; γ k lÀ 1 ; s 2 g � � lÀ 1 ; s 2 e À � lÀ 1 ; m lÀ 1 from the conditional posterior distribution: where: ; s 2 e À � lÀ 1 ; m lÀ 1 from the conditional posterior distribution: where: Step "D" refers to the component of variance related to the singular value for which three hypotheses are specified. Here, only hypothesis "c" will be used to illustrate the sampling process. For the other hypotheses, the process is analogous.
The size of the model is given by the number of multiplicative terms retained to explain GEI (0, 1, 2,. . ., rank(GEI)). For each proposed model, this dimension is directly related to steps C, D, E, and F of the Gibbs sampling process. For the model without the presence of the multiplicative part t = 0, steps C, D, E, and F are not included in the sampling.
RJMCMC algorithm. With the sampling process for the parameters of the candidate models already established, the next step is to test the models, using the decision rule, selecting the most probable. For this, the following steps of the algorithm are performed: i. Inclusion movement 1. Propose the addition of a bilinear term in the model (t + 1); 2. Calculate the probability of accepting a(t, t + 1), the inclusion of the bilinear term; 3. Sampling a random number from a uniform distribution u~U(0,1); 4. If u < min(1, a(t, t + 1)), the proposed state is accepted; otherwise, the model is retained with the initial dimension.
ii. Exclusion movement: 1. Propose the removal of a bilinear term in the model (t − 1); 2. Calculate the probability of accepting a(t, t − 1) the exclusion of the bilinear term; 3. Sampling a random number from a uniform distribution u~(0,1); 4. If u < min(1, a(t, t − 1)), the proposed state is accepted; otherwise, the model remains with the initial dimension.
Repeat the steps specified in i) and ii) during the iterative process.
The model selected to represent the dataset will be the one with the highest frequency of visits, based on the decision rule.

The predictive capacity of the model
To assess the predictive capacity of the models with and without the use of the RJMCMC algorithm, an unbalance process was performed on the data such that within each environment, random losses of genotypes were established (total loss of the genotype) at the level of loss of 10% of the data set. Thus, the original sample was randomly divided into 10 sub-samples (10-fold), keeping all the genotype and environment information in the sub-samples. In each fold, the data were divided into training and validation data.
The predictive capacity of the models was evaluated using the PRESS (Predicted Residual Sum of Squares) criteria and the correlation (Cor) between the predicted (ŷ ijr ) and observed (y ijr ) phenotypic values. This verification was given in scenarios involving the different hypotheses of the variance component related to the singular value.
Cor and PRESS are calculated, respectively, by: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P n i¼1ŷijr À � y ijr � � 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P n i¼1 y ijr À � y ijr where � y ijr is the predicted mean values, � y ijr the mean of the observed values and n denotes the number of data removed for validation.

Diagnosis of Markov chains and a posteriori inference
The total number of iterations to be performed with the Gibbs sampler in each analysis, was determined by a pilot sample according to Raftery and Lewis [37], indicating the number of initial observations to be discarded (burning), and the jumps necessary to collect the observations (tinning) that will compose the sample for inference. These procedures were performed to avoid selecting values from non-stationary strings and/or autocorrelated strings. The convergence of the produced Markov chains, after the correction procedures described above, was monitored using the method of Raftery and Lewis [37] and by the criterion of Heidelberger and Welch [38]. Trace and density plots were used for visual analysis of the stationarity of Markov chains.
The singular vectors, in turn, are estimated from a correction to preserve the orthonormality constraints between the vectors [10]. This is done by orthonormalizing the matrices "U" and "V", whose columns are composed by the average coordinates of the genotypic and environmental singular vectors, respectively.
Linear parameters, singular values, and variance components of the selected model were estimated by posterior means of the simulated distributions. In addition, uncertainty regarding parameter estimation was quantified through the highest posterior density credible region (HPD).
Posterior bivariate credibility regions for genotypic and environmental scores were incorporated into the biplot-AMMI using the method described by Hu and Yang [39] to study stability and adaptability.
For models adjusted by the RJMCMC method, posterior means, HPD regions with 95% credibility, and bivariate credibility regions for genotypic and environmental scores in the biplot were obtained from two perspectives: i. Based on the most likely model, which is the one that considers only the observations contained in the conditional model (the model with the highest frequency of "visits"), which is referred to as conditional response herewith.
ii. Based on the marginal model that considers all MCMC observations of the possible models that contain the estimated parameters present in the most likely model; that is, it considers the estimates of nested models in which the number of bilinear terms is equal to or greater than that of the most likely model. This is referred to in the text as marginal response.
The entire inference process, as well as the convergence tests, were conducted using the boa package [40] and other tools from the R statistical software.

Results
The results obtained by applying the criteria of Heidelberger and Welch [38], and Raftery and Lewis [41] indicated good convergence properties for MCMC chains (Gibbs sampling) for all parameters with a dependence factor (I) always lower than five (I<5). Furthermore, all parameters passed the stationarity test, indicating that convergence was achieved.

Model selection using information criteria
Based on Fig 1, the model with the best fit, according to the AICM criterion, is independent of the assumed prior for the variance components of the singular value, being the one with two retained bilinear terms (AMMI-2). For BAMMIE, the model that best fits the data was also the same (AMMI-2), regardless of the selection criteria used. For BAMMI and BAMMIS, the best fit models were, respectively, AMMI-3 for the BIC; for the AIC, it was AMMI-4 and AMMI-3, for BAMMI and BAMMIS, respectively. . It is noteworthy that this information, by itself, does not correspond to tests to determine the number of bilinear components that should be retained in the model, but is useful for describing the behavior of the data and can still be used as auxiliary criteria in the step of selection. Table 1 presents posterior means of the singular values for the models as a function of the number of retained bilinear terms, as well as the respective sum of squares of the accumulated GEI to each dimension. An accentuated shrinkage effect was observed in the estimates of the singular values of BAMMIE, concerning the other models, from the adjustment of the third multiplicative term. This accentuated shrinkage was also seen in BAMMIS concerning BAMMI, as from the model with a three-axis. S1 Table presents "true" values for the λ k , together with the respective GEI sums of squares, conditioned to each dimension of the model, and were obtained by Singular Value Decomposition (SVD) of the simulated interaction matrix, in which the sum of squares accumulated was 259.6064.
A point to be noted is that no member of the BAMMIE family overestimated the total variance of the GEI (simulated), and its estimates were always smaller concerning BAMMIS and BAMMI. Combining information from Fig 1 and Table 1, it is observed that models selected by the AICM criterion did not present an accumulated sum of squares greater than the true total sum of squares. This fact was also verified for BAMMIE, regardless of the criteria. For the models selected in BAMMIS and BAMMI using AIC and BIC, the sum of squared estimates retrieved was greater than the total simulated sum of squares.
Another finding was that the BAMMIE model selected by all information criteria (AMMI-2) presented a less discrepant accumulated sum concerning the first two terms of the SVD solution (S1 Table) compared to the others. However, considering complete models, BAMMIE was the one that most underestimated the true total sum of squares, BAMMI was the one that most overestimated it, and BAMMIS offered the closest value.
Summaries with a posteriori estimates of model parameters, including credibility regions for genotype effects and biplot representations, with incorporated bivariate credibility regions for genotypic and environmental scores that describe the GEI, referring to the adjustment using only Gibbs sampling, are presented in S1-S6 Figs and S2 Table. Adjustment of the AMMI model using the Reversible Jump method Fig 2 illustrates the frequencies at which each model was visited in the RJMCMC sampling. These frequencies were used as criteria to select the winning model (most likely). According to this criterion, models AMMI3, AMMI2, and AMMI3 respectively would be selected for BAMMI, BAMMIE, and BAMMIS. Note that there are significant differences between the frequency of the winning model concerning the others for each AMMI version. Table 2 indicates that the models presented very close estimates for the mean of the first two singular values. Furthermore, BAMMI and BAMMIS obtained similar mean values

PLOS ONE
concerning λ3. Also, in relation to λ 3 , there is a more expressive difference for BAMMIE, although the HPD intervals overlap with 95% probability.
The sum of squares for the selected models BAMMI, BAMMIE, and BAMMIS were 272.7710, 217.3316, and 264.2804, respectively. Thus, the most probable models for BAMMI and BAMMIS overestimated the GEI total sum of squares, while BAMMIE had a relatively smaller value and was closer to the true sum of squares than that obtained by the AICM criterion (S1 Table). The accentuated shrinkage effect for the posterior mean of λ3 in the BAMMIE model, in relation to the others, as well as in relation to the real value (S1 Table), is remarkable, which indicates that the estimates for the other parameters should be shrunk to zero. This sharp shrinkage effect was not observed in the other selected models.
Posterior means and HPD credibility regions for residual and genotypic variances for the conditional and marginal responses of the models are presented in Table 3. No large discrepancies were observed between the conditional and marginal estimates or between the three models. Even when compared with the estimates of the models selected by the information

PLOS ONE
criteria, it was found that there are no major discrepancies about genotypic and residual variances (S2 Table). In Fig 3, posterior means and HPD intervals at 95% of credibility of the genotypic effects for BAMMI, BAMMIS, and BAMMIE are presented considering conditional (blue color) and marginal (red color) responses. In both situations, the pattern in the classification of main

PLOS ONE
genotypic effects was practically the same between the models. Overlap between ranges indicates similar effects. Regarding yield, the same subgroup of genotypes (more productive) {G1, G12, G4, G10, G13, G11} was selected for the three models. The subgroup {G19, G2, G16, G9, G18, G5, G14} had a significantly lower mean than the overall one (negative effect), not being interesting in terms of productivity. Lastly, there are those genotypes whose HPD regions include the origin and, therefore, the effects do not differ significantly from zero. Differences between the main effect estimates between the models are practically imperceptible graphically. These results are presented in S3-S5 Tables.

Biplot credibility regions for models adjusted by Reversible Jump
In Fig 4, biplots with posterior credibility regions for the genotypic scores of the three models are presented, with two representations for each model: marginal response and conditional Table 3. Posterior means the standard deviation (SD) and the 95% HPD intervals (LL: lower limit, LU: upper limit) for components of variance genotypic (σ g 2 ) and residual (σ e 2 ) for models with conditional and response responses.

Model
Par  response (as explained earlier in the Methods). From the visualization of the biplots, it is clear (for each model) that the two representations (marginal or conditional) did not result in differences in interpretation. However, there is clearly a change in the stability rating when comparing different models. As the responses within each model did not differ, conclusions about existing GEI standards will be made based on the marginal response. Fig 5 presents the biplot with the credibility regions for genotypic and environmental scores (at 95% probability) referent the marginal BAMMI model. Regions that do encompass the origin (0,0) are not plotted for ease of interpretation biplot. The subgroup of genotypes formed by {G2, G4, G15, G17, G18, G19, G20}, which were not represented, were interpreted as stable, as they have no important contribution to the GEI. The same interpretation applies to the environments and, in this sense, E5, which was not represented, was considered stable. The subgroups of genotypes and environments, whose associated bivariate regions do not encompass the origin, effectively contribute to the GEI effect and are often referred to as unstable.
Furthermore, it is possible to identify similar subgroups of genotypes and environments concerning the interaction effect (homogeneous subgroups). Although with some overlaps that make it difficult to observe separable groups, one could identify five homogeneous sub- The biplots for the BAMMIS and BAMMIE models are shown in Figs 6 and 7, respectively. The only difference between BAMMI and BAMMIS is the G3 genotype, which was classified as stable in the BAMMIS analysis. In the BAMMIE biplot, the credibility regions of G3 and G16 also encompass the origin, signaling that these genotypes would not be important for the effect of GEI, which are the main differences in interpretations of BAMMIE concerning BAMMI and BAMMIS.

Evaluation of the predictive capacity of the models
The evaluation of the predictive capacity of the models utilizing the Cor (the higher the value, the better) and PRESS (smaller values indicate a greater predictive capacity) statistics, was performed using the Gibbs and RJMCMC (marginal responses of the models) algorithms. Fig 8 illustrates the behaviors of the three AMMI versions adjusted by the Gibbs sampler, based on mean Cor and PRESS. The results suggest a high predictive capacity for the models under the used imbalance (always above 81%), with emphasis on the BAMMI-1 that has a correlation above 84%. For the PRESS criterion, the models that presented the best results were BAMMI-1, BAMMIS-1, and BAMMIE-2, similar to the results obtained based on the correlation; note that the BAMMI-1 model achieved slightly higher results in both evaluation criteria.
The results of Cor and PRESS for the models fitted by the RJMCMC method are presented in Fig 9. For Cor, a mean behavior above 89% was observed, and individually (for the folds) the results were always above 84% for all models. For the PRESS criterion, the models that presented the best results were BAMMI-4, BAMMIS-4 and BAMMIE-8, and for Cor, the best result was presented by AMMI-8, regardless of the model. A similar pattern was observed between BAMMI and BAMMIS, both of which were better (in prediction) than BAMMIE, which also showed good performance.

Real data
Analyses of actual data were performed using only the RJMCMC. The posterior estimates were obtained considering the marginal and conditional responses of the models, in a similar manner to what was done with the simulated data. In Fig 10, the histograms of the frequency of visits of the models in relation to the number of bilinear terms are presented. Based on the frequencies of visits, AMMI-3 would be selected for the considered three versions of the AMMI model. Similar to in the study of simulated data, differences were observed in the frequency of the most likely model concerning the others. Table 4 presents posterior means of singular values (for the conditional and marginal responses) and credibility regions estimated based on the marginal response. There were no accented differences between point estimates (conditional and marginal). Among the models, the mean estimates were close for the first three singular values, a fact that was also observed for the HPD interval at 95% of credibility.
Posterior means along with HPD regions for the residual and genotypic variance components, for the marginal and conditional responses, are presented in Table 5. The estimates, considering conditional and marginal responses, are practically the same within each model. Between the models, there were also no expressive differences in the estimates of the corresponding parameters.

PLOS ONE
In Fig 11, posterior means and HPD intervals at 95% of credibility are presented, for the effect of genotypes, conditional response (blue color), and marginal response (red color). In both situations, the pattern in the classification of the effects was practically the same within each model. The inferences were also practically the same in both models with the same subgroup of genotypes (more productive) being selected {G9, G10, G38, G8, G6, G2, G1, G3, G5, G7}. , are presented as biplots with 95% bivariate credibility regions for genotypic and environmental scores' (to the AMMI3 model) marginal response. Here too, only the regions   The adaptability of genotypes (or subgroup of genotypes) to environments (or subgroups of environments) can also be suggested from the biplot visualization. This analysis can be performed by observing positions and overlaps between regions of credibility for genotypic and environmental scores in the quadrants referring to the PC1 and PC2 axes. It is noteworthy that these inferences would be almost identical in both biplots.

Points considered for exemplifying the method
In the introduction, we pose two fundamental questions concerning the work of plant breeders. We previously described some differences between the three versions of the AMMI model fitted with MCMC and RJMCMC. We will now highlight differences (or similarities) between the two adjustment methods (MCMC × RJMCMC) for simulated and real data in the analysis of adaptability and phenotypic stability.

Stability analysis for simulated data
Estimates for the singular values of selected models (from information criteria) and the most likely models using RJMCMC were not very outliers. However, as already reported, it was observed that in the assessment of the predictive capacity, the models adjusted by the RJMCMC method performed better than those adjusted by the MCMC method.
The selection and recommendation of genotypes are one of the main objectives of plant breeding programs. The identification of stable genotypes makes it possible to minimize the interaction effect if they have interesting yields. We can compare the main effects of genotypes between models fitted with MCMC and RJMCMC by looking at Fig 3 and the figures that are in the S1-S3 Figs. From the graphical analysis, we found that the inferences do not change even between the different AMMI versions. The best genotypes in terms of the main effect are the same {G1, G12, G4, G10, G13, G11}.

PLOS ONE
The second step in this evaluation would be to observe which of these genotypes would be stable by comparing the versions of the models with MCMC and RJMCMC. Considering MCMC models (S4-S6 Figs), G4 is stable in all the model AMMI versions, while G13 is stable also for BAMMIS and G12 for BAMMIE. Therefore, among the best genotypes (considering main effects), these would have a broad recommendation, respectively, for BAMMI, BAMMIS, and BAMMIE.

PLOS ONE
When analyzing the biplots of the models for the RJMCMC analyses, it is noticed that G4 is also stable in all of them. The only difference is that G13 is also considered stable for BAMMI-E, which differs from the analysis using MCMC. But, in general, the inferences between these two approaches converge.

Adaptability to simulated data
Another way to mitigate the effect of the GEI is to take advantage of its positive effect. For this, it is necessary to evaluate the biplot and interpret positions and overlaps of genotypic and environmental credibility regions.
In the biplots of models fitted with MCMC (S4-S6 Figs), the visual analysis suggests separable subgroups concerning the interaction effect. For the three versions of the model, the GEI

PLOS ONE
standard, in general, remains in the biplot representation. It is observed that for models fitted with RJMCMC the patterns are similar (Figs 5-7). The similarity is also observed when comparing the biplots of equivalent models (MCMC × RJMCMC). This shows us that the analysis of adaptability and stability when considering the different methods, in general, offers the same inferences. Concerning prediction, a slight advantage was seen for the RJMCMC method. In addition, the computational analysis time for the adjustment proved to be much more feasible.

Stability analysis for real data
We exemplified the RJMCMC method with real data and observed, concerning the main effect, that the subgroup {G9, G10, G38, G8, G6, G2, G1, G3, G5, G7} stood out, and even though this inference is independent of the model (BAMMI, BAMIS or BAMMIE). In addition, we saw that G7 was the one that presented the highest value for the main effect since its HPD region does not overlap with any other (Fig 11).
In the biplot analysis (Figs 12, 13, and 14), we also observed that the genotypes G1, G7, 10, G11, G8, and G39 have significant contributions to the GEI the others and were considered stable since their bivariate credibility regions include the origin (0,0). This inference was common to both versions (BAMMI, BAMMIS, and BAMMIE). Thus, among the best genotypes for the main effect, the subgroup {G2, G3, G5, G6, G9, G38}, formed by stable genotypes, has a broad recommendation for all environments tested.
For comparison purposes, fitted the models by the MCMC method and used information criteria for model selection. Fig 15 presents the results of model selections obtained by AIC, BIC, and AICM. As additional criteria, the behaviors of log-likelihood and residual variance are presented. In this sense, models with lower residual variance are preferable, while higher log-likelihood indicates the best model.
The AICM criterion indicated that the BAMMIS and BAMMIE models with three bilinear terms (BAMMIS-3 and BAMMIE-3) would be the best, while for the BAMMI four principal components (BAMMI-4) would be necessary to obtain the optimal model. Using BIC and AIC criteria, we would select different models for each AMMI version. For the AIC criterion, we would have BAMMIE-3, BAMMIS-5, and BAMMI-6; while the BIC criterion pointed out the models: BAMMIE-3, BAMMIS-4, and BAMMI-5.
In Figs 16-18 present posterior means and HPD regions, at 95% credibility, for genotype main effects of the models BAMMIE-3, BAMMIS-3, and BAMMI-4, respectively. To make these inferences, we considered the models selected by the AICM criterion, remembering that for the BAMMIE, the model selected by the criteria used was BAMMIE-3.
As you can see, the inferences regarding the main effect (ranking and credibility regions) do not change from one version to another. When we compare these results with those obtained with the models adjusted with RJMCMC (Fig 11), the inferences are the same, with the subgroup {G9, G10, G38, G8, G6, G2, G1, G3, G5, G7} highlighting and G7 being the best genotype. This behavior was also observed for models adjusted using AIC or BIC (S4 Appendix).
The biplots for the models selected from the AICM method are shown in Figs 19-21. We observe that the positions and overlaps of the credibility regions in the biplots of the three models are practically identical. Furthermore, there are no significant discrepancies when comparing these biplots with those of Figs 12-14, built from the fit using the RJMCMC method. So, the inferences are also the same.

Adaptability in real data
For models fitted with the RJMCMC method, we highlighted the following separable homogeneous subgroups of unstable genotypes {G1}, {G7, 10}, {G11} and {G8, G39}. Genotypes in each subgroup are similar concerning the effect of GEI. We point out that the same subgroups can be identified regardless of the model (BAMMI, BAMMIS, and BAMMIE). Therefore, considering the RJMCMC biplot analysis (Figs 12-14), it is possible to make the following recommendations: G1 genotype for the E8 environment, the subgroup of genotypes {G8, G39} for the subgroup of environments {E1, E2, E3, E4, E5} and subgroup {G7, G10} to environments E7 and E10.
This same result can be inferred from the biplot analysis for the model adjusted by the MCMC method. As already highlighted, the advantages of using the RJMCMC are related to the probabilistic nature of this proposal, such as its feasibility, since the use of information criteria in the selection of models is exhaustive and may become unfeasible. Concerning the analysis of stability and adaptability, we did not find significant differences when considering the different AMMI versions, depending on the a priori distribution adopted for the singular values. Our objective, in this sense, was to show that the RJMCMC can be applied to any Bayesian version of the AMMI model present in the literature.

Discussion
In analyzes using linear-bilinear models, it is important to determine how many bilinear terms should be retained to describe the pattern of interaction. In the frequentist method, this is usually done through parametric tests based on the F distribution or using cross-validation

PLOS ONE
schemes [5,17,19]. The question of the dimensionality of the model naturally arises from the Bayesian perspective, although this subject is still not very widespread in the Bayesian analysis of linear-bilinear models and, in particular, in the Bayesian-AMMI. Most approaches have not focused on the selection of models, but the assessment of stability and adaptability, using credibility regions incorporated in the biplot [11,13,42,43].
The most recent works on model selection in the Bayesian analysis of AMMI and GGE were carried out by Silva et al. [24] and Oliveira et al. [27], respectively. The analyzes presented by these authors can be considered Bayesian versions of the frequentist proposal of Cornelius and Crossa [25]. One of the reasons for using shrinkage estimators is the possibility of dispensing with the use of any criteria to determine the number of bilinear parameters to be retained in the analysis. The use of shrinkage estimators provides better adjustments in terms of prediction. But the method does not necessarily lead to a more parsimonious choice in terms of the model dimension. In this sense, we argue in favor of the need to think about procedures for selecting models, even when specific priors are used for the components of singular values, as is the case of the BAMMIS and BAMMIE models that were discussed in this study. Liu [10] offers a rich discussion on the selection of models for the Bayesian-AMMI, using information criteria BIC and AIC, and Bayes factor to determine the best model.
Despite their popularity, information criteria such as AIC and BIC require an exhaustive evaluation of all possible models. This fact can make applications unfeasible when the model

PLOS ONE
space is large and, therefore, only in limited circumstances are indicated. Furthermore, the Bayes factor, as Liu [10] argues, is sensitive to the choice of priors, and choosing an ideal prior is a difficult task (although this criterion has not been used in this study). It is also noted that different criteria, in general, tend to choose different models (Fig 1). In this sense, it seems more reasonable to opt for the AICM criterion, which, for the conditional inference, proved to be more consistent, for selecting models with the same dimension. In turn, the use of the RJMCMC algorithm, allows both the adjustment and determination of the model dimension simultaneously. It is one of the most popular methods involving transdimensional Markov chains and has been applied in several contexts [29,34,[44][45][46].
In assessing the predictive ability, we consistently observed good properties for the three versions of the AMMI model. As noted, the results of the models fitted by the two methods do not coincide. The results obtained with the RJMCMC were always superior to the Gibbs approach both by the PRESS criterion and by the Correlation. This was also observed for the conditional response of each model (Figs 8 and 9 and S6 Table).
Another point to be highlighted was the use of different priors related to the parameter s 2 l k . In the BAMMI model, as in Crossa et al. [11] and Oliveira et al. [13], a prior can be considered as an approximation of the principle of insufficient reason. The BAMMIS proposed by Silva et al. [24] was obtained by assigning the Jeffreys prior [47]. This prior is based on the principle of invariance. Other suppositions were investigated by Liu [10], where he showed that different suppositions could lead to different results, as verified in this study as well. Despite the wide use of non-informative priors, they may not represent the current state of available information about the problem. In most cases, they are used only as reference distributions, that is, as a default option in the absence of any available information. Thus, among non-informative prior calls, there are some that are more useful and, therefore, frequently used. But depending on the situation, it would be incorrect to say that these are less informative than others. In this sense, the use of a prior based on the principle of maximum entropy [27] provided us with a theoretical justification for performing inference with less informative priors. Using maximum entropy priors can also be useful in situations of incomplete or doubtful information [48][49][50][51].
With simulated data, we observed that, in general, the choice of priors did not affect the pattern of the marginal or conditional response of the model in the biplot. The prior based on the principle of maximum entropy (BAMMIE) provided a more parsimonious choice, but, with few exceptions, it did not change the patterns of the marginal bivariate credibility regions obtained with other priors. This choice is convenient in the example, as the breeder tends to prefer biplot representations to evaluate genotype responses in different environments. The biplot analysis facilitates the work of selection and recommendation of superior cultivars. Although no expressive discrepancies were found, opting for the marginal biplot could be interesting as it produces less Monte Carlo error. However, the conditional biplot determines a greater degree of freedom for the error.

PLOS ONE
From S8 and S9 Tables, it is possible to infer that the models, in general, performed well in correctly identifying stable and unstable genotypes (based on the simulated scenario). This behavior was verified for the fit using only the Gibbs sampler as well as the RJMCMC method. For real data, the most likely model was the one with three retained multiplicative terms (AMMI3) for all AMMI versions (BAMMI, BAMMIS, and BAMMIE). This became more evident for BAMMIE, as indicated by noting the marked difference between the frequency of "visits" (Fig 10). Furthermore, the biplot patterns for the three models are practically

PLOS ONE
indistinguishable. Therefore, in theory, for a graphical evaluation of patterns, BAMMIE could be preferable. The justification is based on the singular value decomposition property that the first axes capture more information. BAMMIE better estimates the first two singular values, reducing the estimates of those related to higher dimensions to zero. Furthermore, the most popular biplot analysis is the AMMI-2 biplot analysis, which is systematically used even when the selected model dimension is greater than two. In addition, the use of maximum entropy priors makes the analysis more flexible, justifying the shrinkage effect without the need of

PLOS ONE
imposing restrictions on the prior to avoid obtaining inappropriate marginal posterior distributions, as presented in Silva et al. [24]. This aspect was also emphasized by Oliveira et al. [27] who reported that the use of the principle of maximum entropy brings greater stability to the algorithm and simplifies the estimation of parameters.
The issues raised here present a good reason to choose the RJMCMC, although the dataset used in our approach is small. In this sense, the computational time of analysis with RJMCMC is much more attractive. As can be seen from S7 Table, the selection of models by information

PLOS ONE
criteria involves procedures whose total computational time (in minutes) is at least four times longer compared to the RJMCMC method. Another justification for using RJMCMC, as mentioned in this text, would be that this method leads to a natural extension of the Markov chain theory.
Other RJMCMC approaches could be used. We chose the one that seemed the simplest to us and didn't bother comparing different methods trying to find the best one. Our objective was to implement RJMCMC directly in the model fit and to show that the method is a more viable alternative than the use of information criteria. However, the comparison of different RJMCMC methods, aiming to find the best one, is very interesting and can lead to greater efficiency in the adjustment. Another option would be to use non-MCMC methods, which could solve problems such as poor mixing and the difficulty of diagnosing convergence. In some situations, there is evidence that non-MCMC methods can be more efficient than MCMC algorithms [52,53]. Fearnhead [54] emphasizes this idea, reviewing three alternatives to MCMC methods: Importance Sampling, Forward-Backward algorithm, and Sequential Monte Carlo (SMC), discussing the efficiency and limitations of these procedures. The use of non-MCMC procedures could be an option to make the Bayesian analysis of multiplicative models more flexible. For this, further studies are needed and will be the subject of future research.